Introduction

This notebook shows how to plot an XRD plot for the two polymorphs of CsCl ($Pm\overline{3}m$ and $Fm\overline{3}m$) by computing the structure factors.



In [1]:

    
# Set up some imports that we will need
import math, cmath
from pymatgen import Lattice, Structure
from pymatgen.util.plotting_utils import get_publication_quality_plot
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from IPython.display import Image, display
import collections

%matplotlib inline

We will now set up some data and variables.



In [2]:

    
# Atomic scattering factors for Cs and Cl obtained from Structure of Materials.
data = {
    "Cs": [55, 6.062, 155.837, 5.986, 19.695, 3.303, 3.335, 1.096, 0.379],
    "Cl": [17, 1.452, 30.935, 2.292, 9.980, 0.787, 2.234, 0.322, 0.323]
}

CuKa = 0.1542 * 10 # Angstrom

Let us now write the actual method to plot the XRD. This method is a simpler replica of the version that Prof Ong has implemented in the Python Materials Genomics (www.pymatgen.org) materals analysis code. It is less efficient, but the steps are more clearly laid out.



In [3]:

    
def plot_xrd(structure):
    """
    Let's define the XRD plot generation in a function that works for any given structure.
    This takes into account the Lorentz polarization factor and the multiplicity factors, 
    but not the Debye-Waller factor. 
    """
    latt = structure.lattice
    
    # Obtain crystallographic reciprocal lattice and points within
    # limiting sphere (within 2/wavelength)
    recip_latt = latt.reciprocal_lattice_crystallographic
    recip_pts = recip_latt.get_points_in_sphere([[0, 0, 0]], [0, 0, 0], 2 / CuKa)
    
    intensities = {} 
    two_thetas = []
    
    for hkl, g_hkl, ind in sorted(recip_pts, key=lambda i: (i[1], -i[0][2], -i[0][1], -i[0][0])):
        if g_hkl != 0:
            theta = math.asin(CuKa / (2 / g_hkl))
            s = g_hkl / 2
            s_2 = s ** 2
            
            # Calculate the structure factor, given by the sum of the atomic scattering factors.
            f_hkl = 0
            for site in structure:
                el = site.specie.symbol
                d = data[el]
                #Atomic scattering factor equation using parameters for Cs and Cl
                fs = d[0] - 41.78214 * s_2 * sum([d[2 * i + 1] * math.exp(-d[2 * i + 2] * s_2) for i in xrange(4)])
                f_hkl += fs * cmath.exp(2j * math.pi * np.dot(hkl, site.frac_coords))
            
            # Intensity is given by modulus squared of structure factor
            i_hkl = (f_hkl * f_hkl.conjugate()).real

            #This corrects for the Lorentz factor.
            lorentz_factor = (1 + math.cos(2 * theta) ** 2) / (math.sin(theta) ** 2 * math.cos(theta))

            two_theta = 2 * theta / math.pi * 180

            #Deal with floating point precision issues.
            ind = np.where(np.abs(np.subtract(two_thetas, two_theta)) <
                           1e-5)
            if len(ind[0]) > 0:
                intensities[two_thetas[ind[0]]][0] += i_hkl * lorentz_factor
            else:
                intensities[two_theta] = [i_hkl * lorentz_factor, tuple(hkl)]
                two_thetas.append(two_theta)

    max_intensity = max([v[0] for v in intensities.values()])
    plt = get_publication_quality_plot(16, 10)
    for k in sorted(intensities.keys()):
        if k < 90: # Let's limit the plot to < 90 degrees
            v = intensities[k]
            i = v[0] / max_intensity * 100
            plt.plot([k, k], [0, i], color='k', linewidth=3, label=str(v[1]))
            plt.annotate(str(v[1]), xy=[k, i], xytext=[k, i], fontsize=16)
    plt.xlabel("2 \\theta (degrees)")
    plt.ylabel("Intensities (scaled)")
    plt.tight_layout()

$\alpha$-CsCl ($Pm\overline{3}m$)

Let's start with the typical $\alpha$ form of CsCl.



In [4]:

    
# Create CsCl structure
a = 4.209 #Angstrom
latt = Lattice.cubic(a)
structure = Structure(latt, ["Cs", "Cl"], [[0, 0, 0], [0.5, 0.5, 0.5]])
plot_xrd(structure)

Compare it with the experimental XRD pattern below.



In [5]:

    
display(Image(filename=('./PDF - alpha CsCl.png')))

$\beta$-CsCl ($Fm\overline{3}m$)

Let's now look at the $\beta$ (high-temperature) form of CsCl.



In [6]:

    
# Create CsCl structure
a = 6.923 #Angstrom
latt = Lattice.cubic(a)
structure = Structure(latt, ["Cs", "Cs", "Cs", "Cs", "Cl", "Cl", "Cl", "Cl"], 
                      [[0, 0, 0], [0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], 
                       [0.5, 0.5, 0.5], [0, 0, 0.5], [0, 0.5, 0], [0.5, 0, 0]])

plot_xrd(structure)

Compare it with the experimental XRD pattern below.



In [7]:

    
display(Image(filename=('./PDF - beta CsCl.png')))



In [7]: